module mod_clm_snicar
  !
  ! Calculate albedo of snow containing impurities
  ! and the evolution of snow effective radius
  !
  use mod_intkinds
  use mod_realkinds
  use mod_mpmessage
  use mod_mppparam
  use mod_dynparam
  use mod_runparams
  use mod_stdio
  use mod_clm_varcon, only : rpi
  use mod_clm_nchelper

  implicit none

  private

  save

  public :: SNICAR_RT ! Snow albedo and vertically-resolved solar absorption
  public :: SnowAge_grain    ! Snow effective grain size evolution
  public :: SnowAge_init     ! Initial read in of snow-aging file
  public :: SnowOptics_init  ! Initial read in of snow-optics file

  ! minimum allowed snow effective radius (also "fresh snow" value) [microns]
  real(rk8), public, parameter :: snw_rds_min = 54.526_rk8
  ! number of aerosol species in snowpack (indices described above) [nbr]
  integer(ik4),  public, parameter :: sno_nbr_aer =   8
  ! parameter to include organic carbon (OC) in snowpack radiative calculations
  logical,  public, parameter :: DO_SNO_OC =    .false.
  ! parameter to include aerosols in snowpack radiative calculations
  logical,  public, parameter :: DO_SNO_AER =   .true.

  ! scavenging factor for hydrophillic BC inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_bcphi = 0.20_rk8
  ! scavenging factor for hydrophobic BC inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_bcpho = 0.03_rk8
  ! scavenging factor for hydrophillic OC inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_ocphi = 0.20_rk8
  ! scavenging factor for hydrophobic OC inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_ocpho = 0.03_rk8
  ! scavenging factor for dust species 1 inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_dst1  = 0.02_rk8
  ! scavenging factor for dust species 2 inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_dst2  = 0.02_rk8
  ! scavenging factor for dust species 3 inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_dst3  = 0.01_rk8
  ! scavenging factor for dust species 4 inclusion in meltwater [frc]
  real(rk8), public, parameter :: scvng_fct_mlt_dst4  = 0.01_rk8

  ! Aerosol species indices:
  !  1= hydrophillic black carbon
  !  2= hydrophobic black carbon
  !  3= hydrophilic organic carbon
  !  4= hydrophobic organic carbon
  !  5= dust species 1
  !  6= dust species 2
  !  7= dust species 3
  !  8= dust species 4
  ! number of spectral bands used in snow model [nbr]
  integer(ik4),  parameter :: numrad_snw  =   5
  ! first band index in near-IR spectrum [idx]
  integer(ik4),  parameter :: nir_bnd_bgn =   2
  ! ending near-IR band index [idx]
  integer(ik4),  parameter :: nir_bnd_end =   5

  ! number of effective radius indices used in Mie lookup table [idx]
  integer(ik4),  parameter :: idx_Mie_snw_mx = 1471
  ! maxiumum temperature index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_T_max      = 11
  ! minimum temperature index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_T_min      = 1
  ! maxiumum temperature gradient index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_Tgrd_max   = 31
  ! minimum temperature gradient index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_Tgrd_min   = 1
  ! maxiumum snow density index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_rhos_max   = 8
  ! minimum snow density index used in aging lookup table [idx]
  integer(ik4),  parameter :: idx_rhos_min   = 1

  ! maximum effective radius defined in Mie lookup table [microns]
  integer(ik4),  parameter :: snw_rds_max_tbl = 1500
  ! minimium effective radius defined in Mie lookup table [microns]
  integer(ik4),  parameter :: snw_rds_min_tbl = 30
  ! maximum allowed snow effective radius [microns]
  real(rk8), parameter :: snw_rds_max     = 1500._rk8
  ! effective radius of re-frozen snow [microns]
  real(rk8), parameter :: snw_rds_refrz   = 1000._rk8

  ! minimum snow mass required for SNICAR RT calculation [kg m-2]
  real(rk8), parameter :: min_snw = 1.0e-30_rk8

  ! constant for liquid water grain growth [m3 s-1], from Brun89
  !real(rk8), parameter :: C1_liq_Brun89 = 1.28e-17_rk8
  ! constant for liquid water grain growth [m3 s-1], from Brun89:
  !   zeroed to accomodate dry snow aging
  real(rk8), parameter :: C1_liq_Brun89 = 0._rk8
  ! constant for liquid water grain growth [m3 s-1], from Brun89:
  !   corrected for LWC in units of percent
  real(rk8), parameter :: C2_liq_Brun89 = 4.22e-13_rk8

  ! time constant for removal of BC in snow on sea-ice [s-1]
  !  (50% mass removal/year)
  real(rk8), parameter :: tim_cns_bc_rmv  = 2.2e-8_rk8
  ! time constant for removal of OC in snow on sea-ice [s-1]
  !  (50% mass removal/year)
  real(rk8), parameter :: tim_cns_oc_rmv  = 2.2e-8_rk8
  ! time constant for removal of dust in snow on sea-ice [s-1]
  !  (50% mass removal/year)
  real(rk8), parameter :: tim_cns_dst_rmv = 2.2e-8_rk8

  ! scaling of the snow aging rate (tuning option):
  ! flag for scaling the snow aging rate by some arbitrary factor
  logical :: flg_snoage_scl    = .false.
  ! arbitrary factor applied to snow aging rate
  real(rk8), parameter :: xdrdt = 1.0_rk8

  ! snow and aerosol Mie parameters:
  ! (arrays declared here, but are set in iniTimeConst)
  ! (idx_Mie_snw_mx is number of snow radii with defined parameters
  ! (i.e. from 30um to 1500um))

  ! direct-beam weighted ice optical properties
  real(rk8) :: ss_alb_snw_drc(idx_Mie_snw_mx,numrad_snw)
  real(rk8) :: asm_prm_snw_drc(idx_Mie_snw_mx,numrad_snw)
  real(rk8) :: ext_cff_mss_snw_drc(idx_Mie_snw_mx,numrad_snw)

  ! diffuse radiation weighted ice optical properties
  real(rk8) :: ss_alb_snw_dfs(idx_Mie_snw_mx,numrad_snw)
  real(rk8) :: asm_prm_snw_dfs(idx_Mie_snw_mx,numrad_snw)
  real(rk8) :: ext_cff_mss_snw_dfs(idx_Mie_snw_mx,numrad_snw)

  ! hydrophiliic BC
  real(rk8) :: ss_alb_bc1(numrad_snw)
  real(rk8) :: asm_prm_bc1(numrad_snw)
  real(rk8) :: ext_cff_mss_bc1(numrad_snw)

  ! hydrophobic BC
  real(rk8) :: ss_alb_bc2(numrad_snw)
  real(rk8) :: asm_prm_bc2(numrad_snw)
  real(rk8) :: ext_cff_mss_bc2(numrad_snw)

  ! hydrophobic OC
  real(rk8) :: ss_alb_oc1(numrad_snw)
  real(rk8) :: asm_prm_oc1(numrad_snw)
  real(rk8) :: ext_cff_mss_oc1(numrad_snw)

  ! hydrophilic OC
  real(rk8) :: ss_alb_oc2(numrad_snw)
  real(rk8) :: asm_prm_oc2(numrad_snw)
  real(rk8) :: ext_cff_mss_oc2(numrad_snw)

  ! dust species 1:
  real(rk8) :: ss_alb_dst1(numrad_snw)
  real(rk8) :: asm_prm_dst1(numrad_snw)
  real(rk8) :: ext_cff_mss_dst1(numrad_snw)

  ! dust species 2:
  real(rk8) :: ss_alb_dst2(numrad_snw)
  real(rk8) :: asm_prm_dst2(numrad_snw)
  real(rk8) :: ext_cff_mss_dst2(numrad_snw)

  ! dust species 3:
  real(rk8) :: ss_alb_dst3(numrad_snw)
  real(rk8) :: asm_prm_dst3(numrad_snw)
  real(rk8) :: ext_cff_mss_dst3(numrad_snw)

  ! dust species 4:
  real(rk8) :: ss_alb_dst4(numrad_snw)
  real(rk8) :: asm_prm_dst4(numrad_snw)
  real(rk8) :: ext_cff_mss_dst4(numrad_snw)

  ! best-fit parameters for snow aging defined over:
  !  11 temperatures from 225 to 273 K
  !  31 temperature gradients from 0 to 300 K/m
  !   8 snow densities from 0 to 350 kg/m3
  ! (arrays declared here, but are set in iniTimeConst)
  ! (idx_rhos_max,idx_Tgrd_max,idx_T_max)
  real(rk8), pointer, contiguous :: snowage_tau(:,:,:)
  ! (idx_rhos_max,idx_Tgrd_max,idx_T_max)
  real(rk8), pointer, contiguous :: snowage_kappa(:,:,:)
  ! idx_rhos_max,idx_Tgrd_max,idx_T_max)
  real(rk8), pointer, contiguous :: snowage_drdt0(:,:,:)

  contains

  !
  ! Determine reflectance of, and vertically-resolved solar absorption in,
  ! snow with impurities.
  !
  ! Original references on physical models of snow reflectance include:
  ! Wiscombe and Warren [1980] and Warren and Wiscombe [1980],
  ! Journal of Atmospheric Sciences, 37,
  !
  ! The multi-layer solution for multiple-scattering used here is from:
  ! Toon et al. [1989], Rapid calculation of radiative heating rates and
  ! photodissociation rates in inhomogeneous multiple scattering atmospheres,
  ! J. Geophys. Res., 94, D13, 16287-16301
  !
  ! The implementation of the SNICAR model in CLM/CSIM is described in:
  ! Flanner, M., C. Zender, J. Randerson, and P. Rasch [2007],
  ! Present-day climate forcing and response from black carbon in snow,
  ! J. Geophys. Res., 112, D11202, doi: 10.1029/2006JD008003
  !
  subroutine SNICAR_RT(flg_snw_ice, lbc, ubc, num_nourbanc, filter_nourbanc, &
                       coszen, flg_slr_in, h2osno_liq, h2osno_ice, snw_rds,  &
                       mss_cnc_aer_in, albsfc, albout, flx_abs)
    use mod_clm_type
    use mod_clm_varpar, only : nlevsno, numrad
    implicit none
    ! flag: =1 when called from CLM, =2 when called from CSIM
    integer(ik4), intent(in) :: flg_snw_ice
    integer(ik4), intent(in) :: lbc, ubc       ! column index bounds [unitless]
    ! number of columns in non-urban filter
    integer(ik4), intent(in) :: num_nourbanc
    ! column filter for non-urban points
    integer(ik4), intent(in) :: filter_nourbanc(ubc-lbc+1)
    ! cosine of solar zenith angle for next time step (col) [unitless]
    real(rk8), intent(in) :: coszen(lbc:ubc)
    ! flag: =1 for direct-beam incident flux, =2 for diffuse incident flux
    integer(ik4), intent(in) :: flg_slr_in
    ! liquid water content (col,lyr) [kg/m2]
    real(rk8), intent(in) :: h2osno_liq(lbc:ubc,-nlevsno+1:0)
    ! ice content (col,lyr) [kg/m2]
    real(rk8), intent(in) :: h2osno_ice(lbc:ubc,-nlevsno+1:0)
    ! snow effective radius (col,lyr) [microns, m^-6]
    integer(ik4),  intent(in) :: snw_rds(lbc:ubc,-nlevsno+1:0)
    ! mass concentration of all aerosol species (col,lyr,aer) [kg/kg]
    real(rk8), intent(in) :: mss_cnc_aer_in(lbc:ubc,-nlevsno+1:0,sno_nbr_aer)
    ! albedo of surface underlying snow (col,bnd) [frc]
    real(rk8), intent(in) :: albsfc(lbc:ubc,numrad)
    ! snow albedo, averaged into 2 bands
    ! (=0 if no sun or no snow) (col,bnd) [frc]
    real(rk8), intent(out) :: albout(lbc:ubc,numrad)
    ! absorbed flux in each layer per unit flux incident on top of snowpack
    ! (col,lyr,bnd) [frc]
    real(rk8), intent(out) :: flx_abs(lbc:ubc,-nlevsno+1:1,numrad)

    ! negative number of snow layers (col) [nbr]
    integer(ik4),  pointer, contiguous :: snl(:)
    ! snow liquid water equivalent (col) [kg/m2]
    real(rk8), pointer, contiguous :: h2osno(:)
    ! corresponding landunit of column (col) [idx] (debugging only)
    integer(ik4),  pointer, contiguous :: clandunit(:)
    ! columns's gridcell index (col) [idx] (debugging only)
    integer(ik4),  pointer, contiguous :: cgridcell(:)
    ! landunit type (lnd) (debugging only)
    integer(ik4),  pointer, contiguous :: ltype(:)
    ! longitude (degrees) (debugging only)
    real(rk8), pointer, contiguous :: londeg(:)
    ! latitude (degrees) (debugging only)
    real(rk8), pointer, contiguous :: latdeg(:)
    ! fraction of ground covered by snow (0 to 1)
    real(rk8), pointer, contiguous :: frac_sno(:)

    !
    ! variables for snow radiative transfer calculations

    ! Local variables representing single-column values of arrays:
    ! negative number of snow layers [nbr]
    integer(ik4) :: snl_lcl
    ! snow effective radius [m^-6]
    integer(ik4) :: snw_rds_lcl(-nlevsno+1:0)
    ! direct beam incident irradiance [W/m2] (set to 1)
    real(rk8):: flx_slrd_lcl(1:numrad_snw)
    ! diffuse incident irradiance [W/m2] (set to 1)
    real(rk8):: flx_slri_lcl(1:numrad_snw)
    ! aerosol mass concentration (lyr,aer_nbr) [kg/kg]
    real(rk8):: mss_cnc_aer_lcl(-nlevsno+1:0,1:sno_nbr_aer)
    ! total column snow mass [kg/m2]
    real(rk8):: h2osno_lcl
    real(rk8):: h2osno_liq_lcl(-nlevsno+1:0)  ! liquid water mass [kg/m2]
    real(rk8):: h2osno_ice_lcl(-nlevsno+1:0)  ! ice mass [kg/m2]
    ! albedo of underlying surface [frc]
    real(rk8):: albsfc_lcl(1:numrad_snw)
    ! single-scatter albedo of ice grains (lyr) [frc]
    real(rk8):: ss_alb_snw_lcl(-nlevsno+1:0)
    ! asymmetry parameter of ice grains (lyr) [frc]
    real(rk8):: asm_prm_snw_lcl(-nlevsno+1:0)
    ! mass extinction coefficient of ice grains (lyr) [m2/kg]
    real(rk8):: ext_cff_mss_snw_lcl(-nlevsno+1:0)
    ! single-scatter albedo of aerosol species (aer_nbr) [frc]
    real(rk8):: ss_alb_aer_lcl(sno_nbr_aer)
    ! asymmetry parameter of aerosol species (aer_nbr) [frc]
    real(rk8):: asm_prm_aer_lcl(sno_nbr_aer)
    ! mass extinction coefficient of aerosol species (aer_nbr) [m2/kg]
    real(rk8):: ext_cff_mss_aer_lcl(sno_nbr_aer)

    ! Other local variables
    ! two-stream approximation type
    !  (1=Eddington, 2=Quadrature, 3=Hemispheric Mean) [nbr]
    integer(ik4) :: aprx_typ
    ! flag to use Delta approximation (Joseph, 1976) (1= use, 0= don't use)
    integer(ik4) :: delta
    ! weights applied to spectral bands, specific to direct and diffuse cases
    ! (bnd) [frc]
    real(rk8):: flx_wgt(1:numrad_snw)

    ! flag: =1 if there is snow, but zero snow layers,
    !     =0 if at least 1 snow layer [flg]
    integer(ik4) :: flg_nosnl
    ! flag: =1 to redo RT calculation if result is unrealistic
    integer(ik4) :: trip
    ! defines conditions for RT redo (explained below)
    integer(ik4) :: flg_dover

    ! temporary snow albedo [frc]
    real(rk8):: albedo
    ! temporary summation variable for NIR weighting
    real(rk8):: flx_sum
    ! snow albedo by band [frc]
    real(rk8):: albout_lcl(numrad_snw)
    ! absorbed flux per unit incident flux at top of snowpack (lyr,bnd) [frc]
    real(rk8):: flx_abs_lcl(-nlevsno+1:1,numrad_snw)

    ! h2o mass (liquid+solid) in snow layer (lyr) [kg/m2]
    real(rk8):: L_snw(-nlevsno+1:0)
    ! snow optical depth (lyr) [unitless]
    real(rk8):: tau_snw(-nlevsno+1:0)
    ! aerosol mass in snow layer (lyr,nbr_aer) [kg/m2]
    real(rk8):: L_aer(-nlevsno+1:0,sno_nbr_aer)
    ! aerosol optical depth (lyr,nbr_aer) [unitless]
    real(rk8):: tau_aer(-nlevsno+1:0,sno_nbr_aer)
    ! cumulative (snow+aerosol) optical depth [unitless]
    real(rk8):: tau_sum
    ! column optical depth from layer bottom to snowpack top (lyr) [unitless]
    real(rk8):: tau_clm(-nlevsno+1:0)
    ! temporary summation of single-scatter albedo of all aerosols [frc]
    real(rk8):: omega_sum
    ! temporary summation of asymmetry parameter of all aerosols [frc]
    real(rk8):: g_sum

    ! weighted optical depth of snow+aerosol layer (lyr) [unitless]
    real(rk8):: tau(-nlevsno+1:0)
    ! weighted single-scatter albedo of snow+aerosol layer (lyr) [frc]
    real(rk8):: omega(-nlevsno+1:0)
    ! weighted asymmetry parameter of snow+aerosol layer (lyr) [frc]
    real(rk8):: g(-nlevsno+1:0)
    ! transformed (i.e. Delta-Eddington) optical depth of
    !  snow+aerosol layer (lyr) [unitless]
    real(rk8):: tau_star(-nlevsno+1:0)
    ! transformed (i.e. Delta-Eddington) SSA of snow+aerosol layer (lyr) [frc]
    real(rk8):: omega_star(-nlevsno+1:0)
    ! transformed (i.e. Delta-Eddington) asymmetry paramater of
    !  snow+aerosol layer (lyr) [frc]
    real(rk8):: g_star(-nlevsno+1:0)

    ! gridcell, column, and landunit indices [idx]
    integer(ik4) :: g_idx, c_idx, l_idx
    ! spectral band index (1 <= bnd_idx <= numrad_snw) [idx]
    integer(ik4) :: bnd_idx
    ! snow effective radius index for retrieving Mie parameters from
    ! lookup table [idx]
    integer(ik4) :: rds_idx
    ! index of bottom snow layer (0) [idx]
    integer(ik4) :: snl_btm
    ! index of top snow layer (-4 to 0) [idx]
    integer(ik4) :: snl_top
    integer(ik4) :: fc       ! column filter index
    integer(ik4) :: i        ! layer index [idx]
    integer(ik4) :: j        ! aerosol number index [idx]
    integer(ik4) :: n        ! tridiagonal matrix index [idx]

    ! direct-beam radiation at bottom of layer interface (lyr) [W/m^2]
    real(rk8):: F_direct(-nlevsno+1:0)
    ! net radiative flux at bottom of layer interface (lyr) [W/m^2]
    real(rk8):: F_net(-nlevsno+1:0)
    ! net absorbed radiative energy (lyr) [W/m^2]
    real(rk8):: F_abs(-nlevsno+1:0)
    ! total absorbed energy in column [W/m^2]
    real(rk8):: F_abs_sum
    ! upward radiative flux at snowpack top [W/m^2]
    real(rk8):: F_sfc_pls
    ! net flux at bottom of snowpack [W/m^2]
    real(rk8):: F_btm_net
    ! net flux at top of snowpack [W/m^2]
    real(rk8):: F_sfc_net
    ! sum of all energy terms; should be 0.0 [W/m^2]
    real(rk8):: energy_sum
    ! direct-beam radiation at bottom of snowpack [W/m^2]
    real(rk8):: F_direct_btm
    ! cosine of solar zenith angle (used locally) [frc]
    real(rk8):: mu_not

    ! counter for number of times through error loop [nbr]
    integer(ik4) :: err_idx
    ! gridcell latitude (debugging only)
    real(rk8):: lat_coord
    ! gridcell longitude (debugging only)
    real(rk8):: lon_coord
    ! underlying surface type (debugging only)
    integer(ik4) :: sfctype

    ! intermediate variables for radiative transfer approximation:
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: gamma1(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: gamma2(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: gamma3(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: gamma4(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: lambda(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: xgamma(-nlevsno+1:0)
    ! two-stream coefficient from Toon et al. (lyr) [unitless]
    real(rk8):: mu_one
    ! tri-diag intermediate variable from Toon et al. (lyr)
    real(rk8):: e1(-nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (lyr)
    real(rk8):: e2(-nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (lyr)
    real(rk8):: e3(-nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (lyr)
    real(rk8):: e4(-nlevsno+1:0)
    ! intermediate variable: upward flux at bottom interface (lyr) [W/m2]
    real(rk8):: C_pls_btm(-nlevsno+1:0)
    ! intermediate variable: downward flux at bottom interface (lyr) [W/m2]
    real(rk8):: C_mns_btm(-nlevsno+1:0)
    ! intermediate variable: upward flux at top interface (lyr) [W/m2]
    real(rk8):: C_pls_top(-nlevsno+1:0)
    ! intermediate variable: downward flux at top interface (lyr) [W/m2]
    real(rk8):: C_mns_top(-nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: A(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: B(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: D(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: E(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: AS(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: DS(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: X(-2*nlevsno+1:0)
    ! tri-diag intermediate variable from Toon et al. (2*lyr)
    real(rk8):: Y(-2*nlevsno+1:0)

    ! Assign local pointers to derived subtypes components (column-level)
    ! (CLM-specific)

    if (flg_snw_ice == 1) then
      snl       => clm3%g%l%c%cps%snl
      h2osno    => clm3%g%l%c%cws%h2osno
      clandunit => clm3%g%l%c%landunit  ! (debug only)
      cgridcell => clm3%g%l%c%gridcell  ! (debug only)
      ltype     => clm3%g%l%itype       ! (debug only)
      londeg    => clm3%g%londeg        ! (debug only)
      latdeg    => clm3%g%latdeg        ! (debug only)
    end if

    frac_sno    => clm3%g%l%c%cps%frac_sno_eff
    clandunit   => clm3%g%l%c%landunit
    ltype       => clm3%g%l%itype

    ! always use Delta approximation for snow
    delta = 1

    ! Loop over all non-urban columns
    ! (when called from CSIM, there is only one column)
    do fc = 1,num_nourbanc
      c_idx = filter_nourbanc(fc)

      ! Zero absorbed radiative fluxes:
      do i=-nlevsno+1,1,1
        flx_abs_lcl(:,:)   = 0._rk8
        flx_abs(c_idx,i,:) = 0._rk8
      end do

      ! set snow/ice mass to be used for RT:
      if (flg_snw_ice == 1) then
        h2osno_lcl = h2osno(c_idx)
      else
        h2osno_lcl = h2osno_ice(c_idx,0)
      end if

      ! Qualifier for computing snow RT:
      !  1) sunlight from atmosphere model
      !  2) minimum amount of snow on ground.
      !     Otherwise, set snow albedo to zero
      if ((coszen(c_idx) > 0._rk8) .and. (h2osno_lcl > min_snw)) then

        ! Set variables specific to CLM
        if (flg_snw_ice == 1) then
          ! Assign local (single-column) variables to global values
          ! If there is snow, but zero snow layers, we must create a
          ! layer locally.
          ! This layer is presumed to have the fresh snow effective radius.
          if (snl(c_idx) > -1) then
            flg_nosnl         =  1
            snl_lcl           =  -1
            h2osno_ice_lcl(0) =  h2osno_lcl
            h2osno_liq_lcl(0) =  0._rk8
            snw_rds_lcl(0)    =  nint(snw_rds_min)
          else
            flg_nosnl         =  0
            snl_lcl           =  snl(c_idx)
            h2osno_liq_lcl(:) =  h2osno_liq(c_idx,:)
            h2osno_ice_lcl(:) =  h2osno_ice(c_idx,:)
            snw_rds_lcl(:)    =  snw_rds(c_idx,:)
          end if

          snl_btm   = 0
          snl_top   = snl_lcl+1

          ! for debugging only
          l_idx     = clandunit(c_idx)
          g_idx     = cgridcell(c_idx)
          sfctype   = ltype(l_idx)
          lat_coord = latdeg(g_idx)
          lon_coord = londeg(g_idx)

          ! Set variables specific to CSIM
        else
          flg_nosnl         = 0
          snl_lcl           = -1
          h2osno_liq_lcl(:) = h2osno_liq(c_idx,:)
          h2osno_ice_lcl(:) = h2osno_ice(c_idx,:)
          snw_rds_lcl(:)    = snw_rds(c_idx,:)
          snl_btm           = 0
          snl_top           = 0
          sfctype           = -1
          lat_coord         = -90
          lon_coord         = 0
        end if

        ! Set local aerosol array
        do j = 1, sno_nbr_aer
          mss_cnc_aer_lcl(:,j) = mss_cnc_aer_in(c_idx,:,j)
        end do

        ! Set spectral underlying surface albedos to their
        ! corresponding VIS or NIR albedos
        albsfc_lcl(1)                       = albsfc(c_idx,1)
        albsfc_lcl(nir_bnd_bgn:nir_bnd_end) = albsfc(c_idx,2)

        ! Error check for snow grain size:
        do i = snl_top, snl_btm, 1
          if ( (snw_rds_lcl(i) < snw_rds_min_tbl) .or. &
               (snw_rds_lcl(i) > snw_rds_max_tbl) ) then
            write (stderr,*) "SNICAR ERROR: snow grain radius of ", &
                    snw_rds_lcl(i), " out of bounds."
            write (stderr,*) "date = ", trim(rcmtimer%str( ))
            write (stderr,*) "flg_snw_ice= ", flg_snw_ice
            write (stderr,*) "column: ", c_idx, " level: ", i, &
                    " snl(c)= ", snl_lcl
            write (stderr,*) "lat= ", lat_coord, " lon= ", lon_coord
            write (stderr,*) "h2osno(c)= ", h2osno_lcl
            call fatal(__FILE__,__LINE__,'clm now stopping')
          end if
        end do

        ! Incident flux weighting parameters
        !  - sum of all VIS bands must equal 1
        !  - sum of all NIR bands must equal 1
        !
        ! Spectral bands (5-band case)
        !  Band 1: 0.3-0.7um (VIS)
        !  Band 2: 0.7-1.0um (NIR)
        !  Band 3: 1.0-1.2um (NIR)
        !  Band 4: 1.2-1.5um (NIR)
        !  Band 5: 1.5-5.0um (NIR)
        !
        ! The following weights are appropriate for surface-incident
        ! flux in a mid-latitude winter atmosphere
        !
        ! 3-band weights
        if ( numrad_snw == 3 ) then
          ! Direct:
          if ( flg_slr_in == 1 ) then
            flx_wgt(1) = 1._rk8
            flx_wgt(2) = 0.66628670195247_rk8
            flx_wgt(3) = 0.33371329804753_rk8
            ! Diffuse:
          else if ( flg_slr_in == 2 ) then
            flx_wgt(1) = 1._rk8
            flx_wgt(2) = 0.77887652162877_rk8
            flx_wgt(3) = 0.22112347837123_rk8
          end if

          ! 5-band weights
        else if (numrad_snw == 5 ) then
          ! Direct:
          if ( flg_slr_in == 1 ) then
            flx_wgt(1) = 1._rk8
            flx_wgt(2) = 0.49352158521175_rk8
            flx_wgt(3) = 0.18099494230665_rk8
            flx_wgt(4) = 0.12094898498813_rk8
            flx_wgt(5) = 0.20453448749347_rk8
            ! Diffuse:
          else if ( flg_slr_in == 2 ) then
            flx_wgt(1) = 1._rk8
            flx_wgt(2) = 0.58581507618433_rk8
            flx_wgt(3) = 0.20156903770812_rk8
            flx_wgt(4) = 0.10917889346386_rk8
            flx_wgt(5) = 0.10343699264369_rk8
          end if
        end if

        ! Loop over snow spectral bands
        do bnd_idx = 1, numrad_snw
          mu_not = coszen(c_idx)  ! must set here, because of error handling
          flg_dover = 1           ! default is to redo
          err_idx   = 0           ! number of times through loop

          do while ( flg_dover > 0 )

            ! DEFAULT APPROXIMATIONS:
            !  VIS:       Delta-Eddington
            !  NIR (all): Delta-Hemispheric Mean
            !  WARNING:   DO NOT USE delta-EDDINGTON FOR NIR DIFFUSE
            !   - this sometimes results in negative albedo
            !
            ! ERROR CONDITIONS:
            ! Conditions which cause "trip", resulting in redo of RT
            ! approximation:
            !   1. negative absorbed flux
            !   2. total absorbed flux greater than incident flux
            !   3. negative albedo
            ! NOTE: These errors have only been encountered in spectral
            ! bands 4 and 5
            !
            ! ERROR HANDLING
            !  1st error (flg_dover=2): switch approximation
            !            (Edd->HM or HM->Edd)
            !  2nd error (flg_dover=3): change zenith angle by 0.02
            !            (this happens about 1 in 10^6 cases)
            !  3rd error (flg_dover=4): switch approximation with new zenith
            !  Subsequent errors: repeatedly change zenith and approximations...

            if ( bnd_idx == 1 ) then
              if ( flg_dover == 2 ) then
                aprx_typ = 3
              else if (flg_dover == 3) then
                aprx_typ = 1
                if ( coszen(c_idx) > 0.5_rk8 ) then
                  mu_not = mu_not - 0.02_rk8
                else
                  mu_not = mu_not + 0.02_rk8
                end if
              else if ( flg_dover == 4 ) then
                aprx_typ = 3
              else
                aprx_typ = 1
              end if
            else
              if ( flg_dover == 2 ) then
                aprx_typ = 1
              else if ( flg_dover == 3 ) then
                aprx_typ = 3
                if ( coszen(c_idx) > 0.5_rk8 ) then
                  mu_not = mu_not - 0.02_rk8
                else
                  mu_not = mu_not + 0.02_rk8
                end if
              elseif ( flg_dover == 4 ) then
                aprx_typ = 1
              else
                aprx_typ = 3
              end if
            end if

            ! Set direct or diffuse incident irradiance to 1
            ! (This has to be within the bnd loop because mu_not is
            ! adjusted in rare cases)
            if ( flg_slr_in == 1 ) then
              ! this corresponds to incident irradiance of 1.0
              flx_slrd_lcl(bnd_idx) = 1._rk8/(mu_not*rpi)
              flx_slri_lcl(bnd_idx) = 0._rk8
            else
              flx_slrd_lcl(bnd_idx) = 0._rk8
              flx_slri_lcl(bnd_idx) = 1._rk8
            end if

            ! Pre-emptive error handling: aerosols can reap havoc on
            ! these absorptive bands.
            ! Since extremely high soot concentrations have a negligible
            ! effect on these bands, zero them.
            if ( ( numrad_snw == 5) .and. &
                 ((bnd_idx == 5).or.(bnd_idx == 4)) ) then
              mss_cnc_aer_lcl(:,:) = 0._rk8
            end if

            if ( (numrad_snw == 3).and.(bnd_idx == 3) ) then
              mss_cnc_aer_lcl(:,:) = 0._rk8
            end if

            ! Define local Mie parameters based on snow grain size
            ! and aerosol species, retrieved from a lookup table.
            if ( flg_slr_in == 1 ) then
              do i = snl_top, snl_btm, 1
                rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                ! snow optical properties (direct radiation)
                ss_alb_snw_lcl(i)      = ss_alb_snw_drc(rds_idx,bnd_idx)
                asm_prm_snw_lcl(i)     = asm_prm_snw_drc(rds_idx,bnd_idx)
                ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_drc(rds_idx,bnd_idx)
              end do
            else if ( flg_slr_in == 2 ) then
              do i = snl_top, snl_btm, 1
                rds_idx = snw_rds_lcl(i) - snw_rds_min_tbl + 1
                ! snow optical properties (diffuse radiation)
                ss_alb_snw_lcl(i)      = ss_alb_snw_dfs(rds_idx,bnd_idx)
                asm_prm_snw_lcl(i)     = asm_prm_snw_dfs(rds_idx,bnd_idx)
                ext_cff_mss_snw_lcl(i) = ext_cff_mss_snw_dfs(rds_idx,bnd_idx)
              end do
            end if

            ! aerosol species 1 optical properties
            ss_alb_aer_lcl(1)        = ss_alb_bc1(bnd_idx)
            asm_prm_aer_lcl(1)       = asm_prm_bc1(bnd_idx)
            ext_cff_mss_aer_lcl(1)   = ext_cff_mss_bc1(bnd_idx)

            ! aerosol species 2 optical properties
            ss_alb_aer_lcl(2)        = ss_alb_bc2(bnd_idx)
            asm_prm_aer_lcl(2)       = asm_prm_bc2(bnd_idx)
            ext_cff_mss_aer_lcl(2)   = ext_cff_mss_bc2(bnd_idx)

            ! aerosol species 3 optical properties
            ss_alb_aer_lcl(3)        = ss_alb_oc1(bnd_idx)
            asm_prm_aer_lcl(3)       = asm_prm_oc1(bnd_idx)
            ext_cff_mss_aer_lcl(3)   = ext_cff_mss_oc1(bnd_idx)

            ! aerosol species 4 optical properties
            ss_alb_aer_lcl(4)        = ss_alb_oc2(bnd_idx)
            asm_prm_aer_lcl(4)       = asm_prm_oc2(bnd_idx)
            ext_cff_mss_aer_lcl(4)   = ext_cff_mss_oc2(bnd_idx)

            ! aerosol species 5 optical properties
            ss_alb_aer_lcl(5)        = ss_alb_dst1(bnd_idx)
            asm_prm_aer_lcl(5)       = asm_prm_dst1(bnd_idx)
            ext_cff_mss_aer_lcl(5)   = ext_cff_mss_dst1(bnd_idx)

            ! aerosol species 6 optical properties
            ss_alb_aer_lcl(6)        = ss_alb_dst2(bnd_idx)
            asm_prm_aer_lcl(6)       = asm_prm_dst2(bnd_idx)
            ext_cff_mss_aer_lcl(6)   = ext_cff_mss_dst2(bnd_idx)

            ! aerosol species 7 optical properties
            ss_alb_aer_lcl(7)        = ss_alb_dst3(bnd_idx)
            asm_prm_aer_lcl(7)       = asm_prm_dst3(bnd_idx)
            ext_cff_mss_aer_lcl(7)   = ext_cff_mss_dst3(bnd_idx)

            ! aerosol species 8 optical properties
            ss_alb_aer_lcl(8)        = ss_alb_dst4(bnd_idx)
            asm_prm_aer_lcl(8)       = asm_prm_dst4(bnd_idx)
            ext_cff_mss_aer_lcl(8)   = ext_cff_mss_dst4(bnd_idx)

            ! 1. snow and aerosol layer column mass (L_snw, L_aer [kg/m^2])
            ! 2. optical Depths (tau_snw, tau_aer)
            ! 3. weighted Mie properties (tau, omega, g)

            ! Weighted Mie parameters of each layer
            do i = snl_top, snl_btm, 1
              L_snw(i)   = h2osno_ice_lcl(i)+h2osno_liq_lcl(i)
              tau_snw(i) = L_snw(i)*ext_cff_mss_snw_lcl(i)

              do j = 1, sno_nbr_aer
                L_aer(i,j)   = L_snw(i)*mss_cnc_aer_lcl(i,j)
                tau_aer(i,j) = L_aer(i,j)*ext_cff_mss_aer_lcl(j)
              end do

              tau_sum   = 0._rk8
              omega_sum = 0._rk8
              g_sum     = 0._rk8

              do j = 1, sno_nbr_aer
                tau_sum    = tau_sum + tau_aer(i,j)
                omega_sum  = omega_sum + (tau_aer(i,j)*ss_alb_aer_lcl(j))
                g_sum = g_sum + &
                        (tau_aer(i,j)*ss_alb_aer_lcl(j)*asm_prm_aer_lcl(j))
              end do

              tau(i) = tau_sum + tau_snw(i)
              omega(i) = (1.0_rk8/tau(i)) * &
                      (omega_sum+(ss_alb_snw_lcl(i)*tau_snw(i)))
              g(i) = (1.0_rk8/(tau(i)*omega(i))) * &
                   (g_sum+ (asm_prm_snw_lcl(i)*ss_alb_snw_lcl(i)*tau_snw(i)))
            end do

            ! delta transformations, if requested
            if ( delta == 1 ) then
              do i = snl_top, snl_btm, 1
                g_star(i) = g(i)/(1.0_rk8+g(i))
                omega_star(i) = ((1.0_rk8-(g(i)**2))*omega(i)) / &
                                (1.0_rk8-(omega(i)*(g(i)**2)))
                tau_star(i)   = (1.0_rk8-(omega(i)*(g(i)**2)))*tau(i)
              end do
            else
              do i = snl_top, snl_btm, 1
                g_star(i)     = g(i)
                omega_star(i) = omega(i)
                tau_star(i)   = tau(i)
              end do
            end if

            ! Total column optical depth:
            ! tau_clm(i) = total optical depth above the bottom of layer i
            tau_clm(snl_top) = 0._rk8
            do i = snl_top + 1, snl_btm, 1
              tau_clm(i) = tau_clm(i-1)+tau_star(i-1)
            end do

            ! Direct radiation at bottom of snowpack:
            if ( (tau_clm(snl_btm)+tau_star(snl_btm))/mu_not > 25.0_rk8 ) then
              F_direct_btm = 0.0_rk8
            else
              F_direct_btm = albsfc_lcl(bnd_idx)*mu_not * &
                  exp(-(tau_clm(snl_btm)+tau_star(snl_btm))/mu_not) * &
                  rpi*flx_slrd_lcl(bnd_idx)
            end if


            ! Intermediates
            ! Gamma values are approximation-specific.

            ! Eddington
            if ( aprx_typ == 1 ) then
              do i = snl_top, snl_btm, 1
                gamma1(i) = (7.0_rk8-(omega_star(i) * &
                        (4.0_rk8+(3.0_rk8*g_star(i)))))/4.0_rk8
                gamma2(i) = -(1.0_rk8-(omega_star(i) * &
                        (4.0_rk8-(3.0_rk8*g_star(i)))))/4.0_rk8
                gamma3(i) = (2.0_rk8-(3.0_rk8*g_star(i)*mu_not))/4.0_rk8
                gamma4(i) = 1.0_rk8-gamma3(i)
                mu_one    = 0.5_rk8
              end do
            ! Quadrature
            else if ( aprx_typ == 2 ) then
              do i = snl_top, snl_btm, 1
                gamma1(i) = (3.0_rk8**0.5_rk8)*(2.0_rk8-(omega_star(i) * &
                        (1.0_rk8+g_star(i))))/2.0_rk8
                gamma2(i) = omega_star(i)*(3.0_rk8**0.5_rk8) * &
                        (1.0_rk8-g_star(i))/2.0_rk8
                gamma3(i) = (1.0_rk8-((3.0_rk8**0.50_rk8)*g_star(i)*mu_not))/2._rk8
                gamma4(i) = 1.0_rk8-gamma3(i)
                mu_one = 1.0_rk8/(3.0_rk8**0.50_rk8)
              end do
            ! Hemispheric Mean
            else if ( aprx_typ == 3 ) then
              do i = snl_top, snl_btm, 1
                gamma1(i) = 2.0_rk8 - (omega_star(i)*(1.0_rk8+g_star(i)))
                gamma2(i) = omega_star(i)*(1.0_rk8-g_star(i))
                gamma3(i) = (1.0_rk8-((3.0_rk8**0.5_rk8)*g_star(i)*mu_not))/2.0_rk8
                gamma4(i) = 1-gamma3(i)
                mu_one    = 0.5_rk8
              end do
            end if

            ! Intermediates for tri-diagonal solution
            do i = snl_top, snl_btm, 1
              lambda(i) = sqrt(abs((gamma1(i)**2) - (gamma2(i)**2)))
              xgamma(i)  = gamma2(i)/(gamma1(i)+lambda(i))
              if ( lambda(i)*tau_star(i) > 21.0_rk8 ) then
                e1(i) = 1.0_rk8
                e2(i) = 1.0_rk8
                e3(i) = xgamma(i)
                e4(i) = xgamma(i)
              else
                e1(i) = 1.0_rk8+(xgamma(i)*exp(-lambda(i)*tau_star(i)))
                e2(i) = 1.0_rk8-(xgamma(i)*exp(-lambda(i)*tau_star(i)))
                e3(i) = xgamma(i) + exp(-lambda(i)*tau_star(i))
                e4(i) = xgamma(i) - exp(-lambda(i)*tau_star(i))
              end if
            end do !enddo over snow layers

            ! Intermediates for tri-diagonal solution
            do i = snl_top, snl_btm, 1
              if ( flg_slr_in == 1 ) then
                if ( (tau_clm(i)+tau_star(i))/mu_not > 21.0_rk8 ) then
                  C_pls_btm(i) = 0.0_rk8
                  C_mns_btm(i) = 0.0_rk8
                else
                  C_pls_btm(i) = (omega_star(i)*rpi*flx_slrd_lcl(bnd_idx)* &
                        exp(-(tau_clm(i)+tau_star(i))/mu_not)*     &
                        (((gamma1(i)-(1.0_rk8/mu_not))*gamma3(i)) +  &
                        (gamma4(i)*gamma2(i))))/((lambda(i)**2) -  &
                        (1.0_rk8/(mu_not**2)))
                  C_mns_btm(i) = (omega_star(i)*rpi*flx_slrd_lcl(bnd_idx)* &
                        exp(-(tau_clm(i)+tau_star(i))/mu_not)*    &
                        (((gamma1(i)+(1.0_rk8/mu_not))*gamma4(i)) + &
                        (gamma2(i)*gamma3(i))))/((lambda(i)**2) - &
                        (1.0_rk8/(mu_not**2)))
                end if
                if ( tau_clm(i)/mu_not > 21.0_rk8 ) then
                  C_pls_top(i) = 0.0_rk8
                  C_mns_top(i) = 0.0_rk8
                else
                  C_pls_top(i) = (omega_star(i)*rpi*flx_slrd_lcl(bnd_idx)* &
                        exp(-tau_clm(i)/mu_not)*(((gamma1(i)-(1.0_rk8/mu_not))* &
                        gamma3(i))+(gamma4(i)*gamma2(i))))/((lambda(i)**2) - &
                        (1.0_rk8/(mu_not**2)))
                  C_mns_top(i) = (omega_star(i)*rpi*flx_slrd_lcl(bnd_idx)* &
                        exp(-tau_clm(i)/mu_not)*(((gamma1(i)+(1.0_rk8/mu_not))* &
                        gamma4(i))+(gamma2(i)*gamma3(i))))/((lambda(i)**2) - &
                        (1.0_rk8/(mu_not**2)))
                end if
              else
                C_pls_btm(i) = 0._rk8
                C_mns_btm(i) = 0._rk8
                C_pls_top(i) = 0._rk8
                C_mns_top(i) = 0._rk8
              end if
            end do

            ! Coefficients for tridiaganol matrix solution
            do i = 2*snl_lcl+1, 0, 1
              ! Boundary values for i=1 and i=2*snl_lcl,
              ! specifics for i=odd and i=even
              if ( i == (2*snl_lcl+1) ) then
                A(i) = 0.0_rk8
                B(i) = e1(snl_top)
                D(i) = -e2(snl_top)
                E(i) = flx_slri_lcl(bnd_idx)-C_mns_top(snl_top)
              else if ( i == 0 ) then
                A(i) = e1(snl_btm)-(albsfc_lcl(bnd_idx)*e3(snl_btm))
                B(i) = e2(snl_btm)-(albsfc_lcl(bnd_idx)*e4(snl_btm))
                D(i) = 0.0_rk8
                E(i) = F_direct_btm-C_pls_btm(snl_btm) + &
                         (albsfc_lcl(bnd_idx)*C_mns_btm(snl_btm))
              else if( mod(i,2) == -1 ) then ! If odd and i>=3 (n=1 for i=3)
                n = floor(i/2.0_rk8)
                A(i) = (e2(n)*e3(n))-(e4(n)*e1(n))
                B(i) = (e1(n)*e1(n+1))-(e3(n)*e3(n+1))
                D(i) = (e3(n)*e4(n+1))-(e1(n)*e2(n+1))
                E(i) = (e3(n)*(C_pls_top(n+1) - &
                     C_pls_btm(n)))+(e1(n)*(C_mns_btm(n)-C_mns_top(n+1)))
              else if ( mod(i,2) == 0 ) then    ! If even and i<=2*snl_lcl
                n = (i/2)
                A(i) = (e2(n+1)*e1(n))-(e3(n)*e4(n+1))
                B(i) = (e2(n)*e2(n+1))-(e4(n)*e4(n+1))
                D(i) = (e1(n+1)*e4(n+1))-(e2(n+1)*e3(n+1))
                E(i) = (e2(n+1)*(C_pls_top(n+1) - &
                     C_pls_btm(n)))+(e4(n+1)*(C_mns_top(n+1)-C_mns_btm(n)))
              end if
            end do

            AS(0) = A(0)/B(0)
            DS(0) = E(0)/B(0)

            do i = -1, (2*snl_lcl+1), -1
              X(i)  = 1.0_rk8/(B(i)-(D(i)*AS(i+1)))
              AS(i) = A(i)*X(i)
              DS(i) = (E(i)-(D(i)*DS(i+1)))*X(i)
            end do

            Y(2*snl_lcl+1) = DS(2*snl_lcl+1)
            do i = (2*snl_lcl+2), 0, 1
              Y(i) = DS(i)-(AS(i)*Y(i-1))
            end do

            ! Downward direct-beam and net flux (F_net) at the base
            ! of each layer:
            do i = snl_top, snl_btm, 1
              if ( (tau_clm(i)+tau_star(i))/mu_not > 25.0_rk8 ) then
                F_direct(i) = 0.0_rk8
              else
                F_direct(i) = mu_not*rpi*flx_slrd_lcl(bnd_idx) * &
                          exp(-(tau_clm(i)+tau_star(i))/mu_not)
              end if
              F_net(i) = (Y(2*i-1)*(e1(i)-e3(i))) + &
                      (Y(2*i)*(e2(i)-e4(i))) + &
                      C_pls_btm(i) - C_mns_btm(i) - F_direct(i)
            end do

            ! Upward flux at snowpack top:
            if ( lambda(snl_top)*tau_star(snl_top) > 25 ) then
              F_sfc_pls = Y(2*snl_lcl+1)*xgamma(snl_top) - &
                 Y(2*snl_lcl+2)*xgamma(snl_top) + C_pls_top(snl_top)
            else
              F_sfc_pls = &
                 (Y(2*snl_lcl+1)*(exp(-lambda(snl_top)*tau_star(snl_top)) + &
                    xgamma(snl_top))) + &
                 (Y(2*snl_lcl+2)*(exp(-lambda(snl_top)*tau_star(snl_top)) - &
                    xgamma(snl_top))) + C_pls_top(snl_top)
            end if

            ! Net flux at bottom = absorbed radiation by underlying surface:
            F_btm_net = -F_net(snl_btm)

            ! Bulk column albedo and surface net flux
            albedo = F_sfc_pls/((mu_not*rpi*flx_slrd_lcl(bnd_idx)) + &
                    flx_slri_lcl(bnd_idx))
            F_sfc_net = F_sfc_pls - ((mu_not*rpi*flx_slrd_lcl(bnd_idx)) + &
                    flx_slri_lcl(bnd_idx))

            trip = 0
            ! Absorbed flux in each layer
            do i = snl_top, snl_btm, 1
              if ( i == snl_top ) then
                F_abs(i) = F_net(i)-F_sfc_net
              else
                F_abs(i) = F_net(i)-F_net(i-1)
              end if
              flx_abs_lcl(i,bnd_idx) = F_abs(i)

              ! ERROR check: negative absorption
              if ( flx_abs_lcl(i,bnd_idx) < -0.00001 ) then
                trip = 1
              end if
            end do

            flx_abs_lcl(1,bnd_idx) = F_btm_net

            if ( flg_nosnl == 1 ) then
              ! If there are no snow layers (but still snow), all absorbed
              ! energy must be in top soil layer
              !flx_abs_lcl(:,bnd_idx) = 0._rk8
              !flx_abs_lcl(1,bnd_idx) = F_abs(0) + F_btm_net
              ! changed on 20070408:
              ! OK to put absorbed energy in the fictitous snow layer
              ! because routine SurfaceRadiation handles the case of no
              ! snow layers. Then, if a snow layer is addded between now
              ! and SurfaceRadiation (called in Hydrology1), absorbed
              ! energy will be properly distributed.
              flx_abs_lcl(0,bnd_idx) = F_abs(0)
              flx_abs_lcl(1,bnd_idx) = F_btm_net
            end if

            !Underflow check (we've already tripped the error condition above)
            do i = snl_top, 1, 1
              if ( flx_abs_lcl(i,bnd_idx) < 0._rk8 ) then
                flx_abs_lcl(i,bnd_idx) = 0._rk8
              end if
            end do

            F_abs_sum = 0._rk8
            do i = snl_top, snl_btm, 1
              F_abs_sum = F_abs_sum + F_abs(i)
            end do

            !ERROR check: absorption greater than incident flux
            ! (should make condition more generic than "1._rk8")
            if ( F_abs_sum > 1._rk8 ) then
              trip = 1
            end if

            !ERROR check:
            if ( (albedo < 0._rk8) .and. (trip == 0) ) then
              trip = 1
            end if

            ! Set conditions for redoing RT calculation
            if ( (trip == 1) .and. (flg_dover == 1) ) then
              flg_dover = 2
            else if ( (trip == 1) .and. (flg_dover == 2) ) then
              flg_dover = 3
            else if ( (trip == 1) .and. (flg_dover == 3) ) then
              flg_dover = 4
            else if( (trip == 1) .and. (flg_dover == 4) .and. &
                     (err_idx < 20) ) then
              flg_dover = 3
              err_idx = err_idx + 1
            else if ( (trip == 1) .and. (flg_dover == 4) .and. &
                      (err_idx >= 20) ) then
              flg_dover = 0
              write(stderr,*) "SNICAR ERROR: FOUND A WORMHOLE."
              write(stderr,*) " STUCK IN INFINITE LOOP!"
              write(stderr,*) " Called from: ", flg_snw_ice
              write(stderr,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(c_idx,0)
              write(stderr,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
              write(stderr,*) "SNICAR STATS: h2osno= ", h2osno_lcl, &
                      " snl= ", snl_lcl
              write(stderr,*) "SNICAR STATS: soot1(0)= ", mss_cnc_aer_lcl(0,1)
              write(stderr,*) "SNICAR STATS: soot2(0)= ", mss_cnc_aer_lcl(0,2)
              write(stderr,*) "SNICAR STATS: dust1(0)= ", mss_cnc_aer_lcl(0,3)
              write(stderr,*) "SNICAR STATS: dust2(0)= ", mss_cnc_aer_lcl(0,4)
              write(stderr,*) "SNICAR STATS: dust3(0)= ", mss_cnc_aer_lcl(0,5)
              write(stderr,*) "SNICAR STATS: dust4(0)= ", mss_cnc_aer_lcl(0,6)
              l_idx     = clandunit(c_idx)
              write(stderr,*) "column index: ", c_idx
              write(stderr,*) "landunit type", ltype(l_idx)
              write(stderr,*) "frac_sno: ", frac_sno(c_idx)
              call fatal(__FILE__,__LINE__,'clm now stopping')
            else
              flg_dover = 0
            end if
          end do !enddo while (flg_dover > 0)

          ! Energy conservation check:
          ! Incident direct+diffuse radiation equals
          ! (absorbed+bulk_transmitted+bulk_reflected)
          energy_sum = (mu_not*rpi*flx_slrd_lcl(bnd_idx)) + &
                  flx_slri_lcl(bnd_idx) - (F_abs_sum + F_btm_net + F_sfc_pls)
          if ( abs(energy_sum) > 0.00001_rk8 ) then
            write (stderr,"(a,e14.7,a,i6,a,i6)") &
               "SNICAR ERROR: Energy conservation error of : ", energy_sum, &
               " at : ", trim(rcmtimer%str( )), " at column: ", c_idx
            call fatal(__FILE__,__LINE__,'clm now stopping')
          end if

          albout_lcl(bnd_idx) = albedo

          ! Check that albedo is less than 1
          if ( albout_lcl(bnd_idx) > 1.0 ) then
            write (stderr,*) &
                 "SNICAR ERROR: Albedo > 1.0 at c: ", c_idx, " at ", &
                 trim(rcmtimer%str( ))
            write (stderr,*) "SNICAR STATS: bnd_idx= ",bnd_idx
            write (stderr,*) "SNICAR STATS: albout_lcl(bnd)= ", &
                    albout_lcl(bnd_idx), " albsfc_lcl(bnd_idx)= ", &
                    albsfc_lcl(bnd_idx)
            write (stderr,*) "SNICAR STATS: landtype= ", sfctype
            write (stderr,*) "SNICAR STATS: h2osno= ", h2osno_lcl, &
                    " snl= ", snl_lcl
            write (stderr,*) "SNICAR STATS: coszen= ", coszen(c_idx), &
                    " flg_slr= ", flg_slr_in
            write (stderr,*) "SNICAR STATS: soot(-4)= ", mss_cnc_aer_lcl(-4,1)
            write (stderr,*) "SNICAR STATS: soot(-3)= ", mss_cnc_aer_lcl(-3,1)
            write (stderr,*) "SNICAR STATS: soot(-2)= ", mss_cnc_aer_lcl(-2,1)
            write (stderr,*) "SNICAR STATS: soot(-1)= ", mss_cnc_aer_lcl(-1,1)
            write (stderr,*) "SNICAR STATS: soot(0)= ", mss_cnc_aer_lcl(0,1)
            write (stderr,*) "SNICAR STATS: L_snw(-4)= ", L_snw(-4)
            write (stderr,*) "SNICAR STATS: L_snw(-3)= ", L_snw(-3)
            write (stderr,*) "SNICAR STATS: L_snw(-2)= ", L_snw(-2)
            write (stderr,*) "SNICAR STATS: L_snw(-1)= ", L_snw(-1)
            write (stderr,*) "SNICAR STATS: L_snw(0)= ", L_snw(0)
            write (stderr,*) "SNICAR STATS: snw_rds(-4)= ", snw_rds(c_idx,-4)
            write (stderr,*) "SNICAR STATS: snw_rds(-3)= ", snw_rds(c_idx,-3)
            write (stderr,*) "SNICAR STATS: snw_rds(-2)= ", snw_rds(c_idx,-2)
            write (stderr,*) "SNICAR STATS: snw_rds(-1)= ", snw_rds(c_idx,-1)
            write (stderr,*) "SNICAR STATS: snw_rds(0)= ", snw_rds(c_idx,0)
            call fatal(__FILE__,__LINE__,'clm now stopping')
          end if
        end do   ! loop over wvl bands

        ! Weight output NIR albedo appropriately
        albout(c_idx,1) = albout_lcl(1)
        flx_sum = 0._rk8
        do bnd_idx= nir_bnd_bgn,nir_bnd_end
          flx_sum = flx_sum + flx_wgt(bnd_idx)*albout_lcl(bnd_idx)
        end do
        albout(c_idx,2) = flx_sum / sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))

        ! Weight output NIR absorbed layer fluxes (flx_abs) appropriately
        flx_abs(c_idx,:,1) = flx_abs_lcl(:,1)
        do i = snl_top, 1, 1
          flx_sum = 0._rk8
          do bnd_idx= nir_bnd_bgn,nir_bnd_end
            flx_sum = flx_sum + flx_wgt(bnd_idx)*flx_abs_lcl(i,bnd_idx)
          end do
          flx_abs(c_idx,i,2) = flx_sum / &
                    sum(flx_wgt(nir_bnd_bgn:nir_bnd_end))
        end do
        ! If snow < minimum_snow, but > 0, and there is sun,
        !  set albedo to underlying surface albedo
      else if ( (coszen(c_idx) > 0._rk8) .and. &
                (h2osno_lcl < min_snw) .and. (h2osno_lcl > 0._rk8) ) then
        albout(c_idx,1) = albsfc(c_idx,1)
        albout(c_idx,2) = albsfc(c_idx,2)
      else ! There is either zero snow, or no sun
        albout(c_idx,1) = 0._rk8
        albout(c_idx,2) = 0._rk8
      end if    ! if column has snow and coszen > 0
    end do    ! loop over all columns
  end subroutine SNICAR_RT
  !
  ! Updates the snow effective grain size (radius).
  ! Contributions to grain size evolution are from:
  !   1. vapor redistribution (dry snow)
  !   2. liquid water redistribution (wet snow)
  !   3. re-freezing of liquid water
  !
  ! Vapor redistribution: Method is to retrieve 3 best-bit parameters that
  ! depend on snow temperature, temperature gradient, and density,
  ! that are derived from the microphysical model described in:
  ! Flanner and Zender (2006), Linking snowpack microphysics and albedo
  ! evolution, J. Geophys. Res., 111, D12208, doi:10.1029/2005JD006834.
  ! The parametric equation has the form:
  ! dr/dt = drdt_0*(tau/(dr_fresh+tau))^(1/kappa), where:
  !   r is the effective radius,
  !   tau and kappa are best-fit parameters,
  !   drdt_0 is the initial rate of change of effective radius, and
  !   dr_fresh is the difference between the current and fresh snow states
  !  (r_current - r_fresh).
  !
  ! Liquid water redistribution: Apply the grain growth function from:
  !   Brun, E. (1989), Investigation of wet-snow metamorphism in respect of
  !   liquid-water content, Annals of Glaciology, 13, 22-26.
  !   There are two parameters that describe the grain growth rate as
  !   a function of snow liquid water content (LWC). The "LWC=0" parameter
  !   is zeroed here because we are accounting for dry snowing with a
  !   different representation
  !
  ! Re-freezing of liquid water: Assume that re-frozen liquid water
  !   into an arbitrarily large effective grain size (snw_rds_refrz).
  !   The phenomenon is observed (Grenfell), but so far unquantified,
  !   as far as I am aware.
  !
  subroutine SnowAge_grain(lbc, ubc, num_snowc, filter_snowc, &
                  num_nosnowc, filter_nosnowc)
    use mod_clm_type
    use mod_clm_varpar      , only : nlevsno
    use mod_clm_varcon      , only : spval
    implicit none
    integer(ik4), intent(in) :: lbc, ubc  ! column bounds
    ! number of column snow points in column filter
    integer(ik4), intent(in) :: num_snowc
    ! column filter for snow points
    integer(ik4), intent(in) :: filter_snowc(ubc-lbc+1)
    ! number of column non-snow points in column filter
    integer(ik4), intent(in) :: num_nosnowc
    ! column filter for non-snow points
    integer(ik4), intent(in) :: filter_nosnowc(ubc-lbc+1)

    ! soil and snow temperature (col,lyr) [K]
    real(rk8), pointer, contiguous :: t_soisno(:,:)
    ! negative number of snow layers (col) [nbr]
    integer(ik4),  pointer, contiguous :: snl(:)
    ! ground temperature (col) [K]
    real(rk8), pointer, contiguous :: t_grnd(:)
    ! layer thickness (col,lyr) [m]
    real(rk8), pointer, contiguous :: dz(:,:)
    ! snow water (col) [mm H2O]
    real(rk8), pointer, contiguous :: h2osno(:)
    ! effective grain radius (col,lyr) [microns, m-6]
    real(rk8), pointer, contiguous :: snw_rds(:,:)
    ! effective grain radius, top layer (col) [microns, m-6]
    real(rk8), pointer, contiguous :: snw_rds_top(:)
    ! liquid water fraction (mass) in top snow layer (col) [frc]
    real(rk8), pointer, contiguous :: sno_liq_top(:)
    ! liquid water content (col,lyr) [kg m-2]
    real(rk8), pointer, contiguous :: h2osoi_liq(:,:)
    ! ice content (col,lyr) [kg m-2]
    real(rk8), pointer, contiguous :: h2osoi_ice(:,:)
    ! snow temperature in top layer (col) [K]
    real(rk8), pointer, contiguous :: snot_top(:)
    ! temperature gradient in top layer (col) [K m-1]
    real(rk8), pointer, contiguous :: dTdz_top(:)
    ! snow on ground after interception (col) [kg m-2 s-1]
    real(rk8), pointer, contiguous :: qflx_snow_grnd_col(:)
    ! excess precipitation due to snow capping [kg m-2 s-1]
    real(rk8), pointer, contiguous :: qflx_snwcp_ice(:)
    ! snow freezing rate (col,lyr) [kg m-2 s-1]
    real(rk8), pointer, contiguous :: qflx_snofrz_lyr(:,:)
    logical, pointer, contiguous :: do_capsnow(:)  ! true => do snow capping
    ! fraction of ground covered by snow (0 to 1)
    real(rk8), pointer, contiguous :: frac_sno(:)

    integer(ik4) :: snl_top         ! top snow layer index [idx]
    integer(ik4) :: snl_btm         ! bottom snow layer index [idx]
    integer(ik4) :: i               ! layer index [idx]
    integer(ik4) :: c_idx           ! column index [idx]
    integer(ik4) :: fc              ! snow column filter index [idx]
    ! snow aging lookup table temperature index [idx]
    integer(ik4) :: T_idx
    ! snow aging lookup table temperature gradient index [idx]
    integer(ik4) :: Tgrd_idx
    ! snow aging lookup table snow density index [idx]
    integer(ik4) :: rhos_idx
    ! temperature at upper layer boundary [K]
    real(rk8) :: t_snotop
    ! temperature at lower layer boundary [K]
    real(rk8) :: t_snobtm
    ! snow temperature gradient (col,lyr) [K m-1]
    real(rk8) :: dTdz(lbc:ubc,-nlevsno:0)
    ! snow aging parameter retrieved from lookup table [hour]
    real(rk8) :: bst_tau
    ! snow aging parameter retrieved from lookup table [unitless]
    real(rk8) :: bst_kappa
    ! snow aging parameter retrieved from lookup table [um hr-1]
    real(rk8) :: bst_drdt0
    ! incremental change in snow effective radius [um]
    real(rk8) :: dr
    ! incremental change in snow effective radius from wet growth [um]
    real(rk8) :: dr_wet
    ! difference between fresh snow r_e and current r_e [um]
    real(rk8) :: dr_fresh
    real(rk8) :: newsnow      ! fresh snowfall [kg m-2]
    real(rk8) :: refrzsnow    ! re-frozen snow [kg m-2]
    real(rk8) :: frc_newsnow  ! fraction of layer mass that is new snow [frc]
    ! fraction of layer mass that is old snow [frc]
    real(rk8) :: frc_oldsnow
    ! fraction of layer mass that is re-frozen snow [frc]
    real(rk8) :: frc_refrz
    ! fraction of layer mass that is liquid water[frc]
    real(rk8) :: frc_liq
    real(rk8) :: rhos  ! snow density [kg m-3]
    ! liquid + solid H2O in snow layer [kg m-2]
    real(rk8) :: h2osno_lyr
    ! column average layer thickness [m]
    real(rk8) :: cdz(-nlevsno+1:0)

    ! Assign local pointers to derived subtypes components (column-level)
    t_soisno           => clm3%g%l%c%ces%t_soisno
    snl                => clm3%g%l%c%cps%snl
    t_grnd             => clm3%g%l%c%ces%t_grnd
    dz                 => clm3%g%l%c%cps%dz
    h2osno             => clm3%g%l%c%cws%h2osno
    snw_rds            => clm3%g%l%c%cps%snw_rds
    h2osoi_liq         => clm3%g%l%c%cws%h2osoi_liq
    h2osoi_ice         => clm3%g%l%c%cws%h2osoi_ice
    snot_top           => clm3%g%l%c%cps%snot_top
    dTdz_top           => clm3%g%l%c%cps%dTdz_top
    snw_rds_top        => clm3%g%l%c%cps%snw_rds_top
    sno_liq_top        => clm3%g%l%c%cps%sno_liq_top
    qflx_snow_grnd_col => clm3%g%l%c%cwf%pwf_a%qflx_snow_grnd
    qflx_snwcp_ice     => clm3%g%l%c%cwf%pwf_a%qflx_snwcp_ice
    qflx_snofrz_lyr    => clm3%g%l%c%cwf%qflx_snofrz_lyr
    do_capsnow         => clm3%g%l%c%cps%do_capsnow
    frac_sno           => clm3%g%l%c%cps%frac_sno_eff

    ! loop over columns that have at least one snow layer
#ifdef STDPAR
    do concurrent ( fc = 1:num_snowc ) local(cdz)
#else
    !$acc parallel loop collapse(1) gang vector private(cdz)
    do fc = 1, num_snowc
#endif
      c_idx = filter_snowc(fc)

      snl_btm = 0
      snl_top = snl(c_idx) + 1
      !$acc loop seq
      do i = snl_top, snl_btm
        cdz(i) = frac_sno(c_idx)*dz(c_idx,i)
      end do

      ! loop over snow layers
      !$acc loop seq
      do i = snl_top, snl_btm, 1
        !
        !**********  1. DRY SNOW AGING  ***********
        !
        h2osno_lyr = h2osoi_liq(c_idx,i) + h2osoi_ice(c_idx,i)

        ! temperature gradient
        if ( i == snl_top ) then
          ! top layer
          t_snotop = t_soisno(c_idx,snl_top)
          t_snobtm = (t_soisno(c_idx,i+1)*dz(c_idx,i) &
                   + t_soisno(c_idx,i)*dz(c_idx,i+1)) &
                  / (dz(c_idx,i)+dz(c_idx,i+1))
        else
          t_snotop = (t_soisno(c_idx,i-1)*dz(c_idx,i) &
                  + t_soisno(c_idx,i)*dz(c_idx,i-1)) &
                  / (dz(c_idx,i)+dz(c_idx,i-1))
          t_snobtm = (t_soisno(c_idx,i+1)*dz(c_idx,i) &
                  + t_soisno(c_idx,i)*dz(c_idx,i+1)) &
                  / (dz(c_idx,i)+dz(c_idx,i+1))
        end if

        dTdz(c_idx,i) = abs((t_snotop - t_snobtm) / cdz(i))

        ! snow density
        rhos = (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i)) / cdz(i)

        ! make sure rhos doesn't drop below 50 (see rhos_idx below)
        rhos = max(50._rk8,rhos)

        ! best-fit table indecies

        T_idx    = nint((t_soisno(c_idx,i)-223.0_rk8) / 5.0_rk8) + 1
        Tgrd_idx = nint(dTdz(c_idx,i) / 10.0_rk8) + 1
        rhos_idx = nint((rhos-50.0_rk8) / 50.0_rk8) + 1

        ! boundary check:
        if (T_idx < idx_T_min) then
          T_idx = idx_T_min
        end if
        if (T_idx > idx_T_max) then
          T_idx = idx_T_max
        end if
        if (Tgrd_idx < idx_Tgrd_min) then
          Tgrd_idx = idx_Tgrd_min
        end if
        if (Tgrd_idx > idx_Tgrd_max) then
          Tgrd_idx = idx_Tgrd_max
        end if
        if (rhos_idx < idx_rhos_min) then
          rhos_idx = idx_rhos_min
        end if
        if (rhos_idx > idx_rhos_max) then
          rhos_idx = idx_rhos_max
        end if

        ! best-fit parameters
        bst_tau   = snowage_tau(rhos_idx,Tgrd_idx,T_idx)
        bst_kappa = snowage_kappa(rhos_idx,Tgrd_idx,T_idx)
        bst_drdt0 = snowage_drdt0(rhos_idx,Tgrd_idx,T_idx)

        ! change in snow effective radius, using best-fit parameters
        dr_fresh = snw_rds(c_idx,i)-snw_rds_min
        dr = (bst_drdt0*(bst_tau/(dr_fresh+bst_tau))** &
                (1.0_rk8/bst_kappa)) * (dtsrf/3600.0_rk8)

        !
        !**********  2. WET SNOW AGING  ***********
        !
        ! We are assuming wet and dry evolution occur simultaneously, and
        ! the contributions from both can be summed.
        ! This is justified by setting the linear offset constant
        ! C1_liq_Brun89 to zero [Brun, 1989]

        ! liquid water faction
        frc_liq = min(0.1_rk8, (h2osoi_liq(c_idx,i) / &
                (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i))))

        !dr_wet = 1E6_rk8*(dtsrf*(C1_liq_Brun89 + &
        !  C2_liq_Brun89*(frc_liq**(3))) / (4*rpi*(snw_rds(c_idx,i)/1E6)**(2)))
        !simplified, units of microns:
        dr_wet = 1.0e18_rk8*(dtsrf*(C2_liq_Brun89*(frc_liq**(3))) / &
                  (4.0_rk8*rpi*snw_rds(c_idx,i)**(2)))
        dr = dr + dr_wet

        !
        !**********  3. SNOWAGE SCALING (TURNED OFF BY DEFAULT)  *************
        !
        ! Multiply rate of change of effective radius by some constant, xdrdt
        if ( flg_snoage_scl ) then
          dr = dr*xdrdt
        end if

        !
        !**********  4. INCREMENT EFFECTIVE RADIUS, ACCOUNTING FOR:  ********
        !               DRY AGING
        !               WET AGING
        !               FRESH SNOW
        !               RE-FREEZING
        !
        ! new snowfall [kg/m2]
        if ( do_capsnow(c_idx) ) then
          newsnow = max(0._rk8, (qflx_snwcp_ice(c_idx)*dtsrf))
        else
          newsnow = max(0._rk8, (qflx_snow_grnd_col(c_idx)*dtsrf))
        end if

        ! snow that has re-frozen [kg/m2]
        refrzsnow = max(0._rk8, (qflx_snofrz_lyr(c_idx,i)*dtsrf))
        ! fraction of layer mass that is re-frozen
        frc_refrz = refrzsnow / h2osno_lyr
        ! fraction of layer mass that is new snow
        if ( i == snl_top ) then
          frc_newsnow = newsnow / h2osno_lyr
        else
          frc_newsnow = 0._rk8
        end if

        if ((frc_refrz + frc_newsnow) > 1._rk8) then
          frc_refrz = frc_refrz / (frc_refrz + frc_newsnow)
          frc_newsnow = 1._rk8 - frc_refrz
          frc_oldsnow = 0._rk8
        else
          frc_oldsnow = 1._rk8 - frc_refrz - frc_newsnow
        end if

        ! mass-weighted mean of fresh snow, old snow, and re-frozen
        ! snow effective radius
        snw_rds(c_idx,i) = (snw_rds(c_idx,i)+dr)*frc_oldsnow + &
                snw_rds_min*frc_newsnow + snw_rds_refrz*frc_refrz
        !
        !**********  5. CHECK BOUNDARIES   ***********
        !
        ! boundary check
        if ( snw_rds(c_idx,i) < snw_rds_min ) then
          snw_rds(c_idx,i) = snw_rds_min
        end if

        if ( snw_rds(c_idx,i) > snw_rds_max ) then
          snw_rds(c_idx,i) = snw_rds_max
        end if

        ! set top layer variables for history files
        if ( i == snl_top ) then
          snot_top(c_idx)    = t_soisno(c_idx,i)
          dTdz_top(c_idx)    = dTdz(c_idx,i)
          snw_rds_top(c_idx) = snw_rds(c_idx,i)
          sno_liq_top(c_idx) = h2osoi_liq(c_idx,i) / &
                  (h2osoi_liq(c_idx,i)+h2osoi_ice(c_idx,i))
        end if
      end do
    end do

    ! Special case: snow on ground, but not enough to have defined a snow layer:
    !   set snw_rds to fresh snow grain size:
    do concurrent ( fc = 1:num_nosnowc )
      c_idx = filter_nosnowc(fc)
      if ( h2osno(c_idx) > 0._rk8 ) then
        snw_rds(c_idx,0) = snw_rds_min
      end if
    end do
  end subroutine SnowAge_grain

  subroutine SnowOptics_init( )
    use mod_clm_varctl     , only : fsnowoptics
    type(clm_filetype) :: ncid                        ! netCDF file id
    character(len= 32) :: subname = 'SnowOptics_init' ! subroutine name
    !
    ! Open optics file:
    if (myid == italk) then
      write(stdout,*) 'Attempting to read snow optical properties .....'
    end if
    call clm_openfile(fsnowoptics,ncid)
    if (myid == italk) then
      write(stdout,*) trim(subname),' : ',trim(fsnowoptics)
    end if

    ! direct-beam snow Mie parameters:
    call clm_readvar(ncid,'ss_alb_ice_drc',ss_alb_snw_drc)
    call clm_readvar(ncid,'asm_prm_ice_drc',asm_prm_snw_drc)
    call clm_readvar(ncid,'ext_cff_mss_ice_drc',ext_cff_mss_snw_drc)

    ! diffuse snow Mie parameters
    call clm_readvar(ncid,'ss_alb_ice_dfs',ss_alb_snw_dfs)
    call clm_readvar(ncid,'asm_prm_ice_dfs',asm_prm_snw_dfs)
    call clm_readvar(ncid,'ext_cff_mss_ice_dfs',ext_cff_mss_snw_dfs)

    ! BC species 1 Mie parameters
    call clm_readvar(ncid,'ss_alb_bcphil',ss_alb_bc1)
    call clm_readvar(ncid,'asm_prm_bcphil',asm_prm_bc1)
    call clm_readvar(ncid,'ext_cff_mss_bcphil',ext_cff_mss_bc1)

    ! BC species 2 Mie parameters
    call clm_readvar(ncid,'ss_alb_bcphob',ss_alb_bc2)
    call clm_readvar(ncid,'asm_prm_bcphob',asm_prm_bc2)
    call clm_readvar(ncid,'ext_cff_mss_bcphob',ext_cff_mss_bc2)

    ! OC species 1 Mie parameters
    call clm_readvar(ncid,'ss_alb_ocphil',ss_alb_oc1)
    call clm_readvar(ncid,'asm_prm_ocphil',asm_prm_oc1)
    call clm_readvar(ncid,'ext_cff_mss_ocphil',ext_cff_mss_oc1)

    ! OC species 2 Mie parameters
    call clm_readvar(ncid,'ss_alb_ocphob',ss_alb_oc2)
    call clm_readvar(ncid,'asm_prm_ocphob',asm_prm_oc2)
    call clm_readvar(ncid,'ext_cff_mss_ocphob',ext_cff_mss_oc2)

    ! dust species 1 Mie parameters
    call clm_readvar(ncid,'ss_alb_dust01',ss_alb_dst1)
    call clm_readvar(ncid,'asm_prm_dust01',asm_prm_dst1)
    call clm_readvar(ncid,'ext_cff_mss_dust01',ext_cff_mss_dst1)

    ! dust species 2 Mie parameters
    call clm_readvar(ncid,'ss_alb_dust02',ss_alb_dst2)
    call clm_readvar(ncid,'asm_prm_dust02',asm_prm_dst2)
    call clm_readvar(ncid,'ext_cff_mss_dust02',ext_cff_mss_dst2)

    ! dust species 3 Mie parameters
    call clm_readvar(ncid,'ss_alb_dust03',ss_alb_dst3)
    call clm_readvar(ncid,'asm_prm_dust03',asm_prm_dst3)
    call clm_readvar(ncid,'ext_cff_mss_dust03',ext_cff_mss_dst3)

    ! dust species 4 Mie parameters
    call clm_readvar(ncid,'ss_alb_dust04',ss_alb_dst4)
    call clm_readvar(ncid,'asm_prm_dust04',asm_prm_dst4)
    call clm_readvar(ncid,'ext_cff_mss_dust04',ext_cff_mss_dst4)

    call clm_closefile(ncid)
    if (myid == italk) then
      write(stdout,*) 'Successfully read snow optical properties'
      ! print some diagnostics:
      if ( debug_level > 3) then
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
              &direct-beam ice, rds=100um: ', &
          ss_alb_snw_drc(71,1), ss_alb_snw_drc(71,2), ss_alb_snw_drc(71,3), &
          ss_alb_snw_drc(71,4), ss_alb_snw_drc(71,5)
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &diffuse ice, rds=100um: ', &
          ss_alb_snw_dfs(71,1), ss_alb_snw_dfs(71,2), ss_alb_snw_dfs(71,3), &
          ss_alb_snw_dfs(71,4), ss_alb_snw_dfs(71,5)
        if (DO_SNO_OC) then
          write (stdout,*) 'SNICAR: Including OC aerosols from snow &
                  &radiative transfer calculations'
        else
          write (stdout,*) 'SNICAR: Excluding OC aerosols from snow &
                  &radiative transfer calculations'
        end if
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &hydrophillic BC: ', &
          ss_alb_bc1(1), ss_alb_bc1(2), ss_alb_bc1(3), ss_alb_bc1(4), &
          ss_alb_bc1(5)
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &hydrophobic BC: ', &
          ss_alb_bc2(1), ss_alb_bc2(2), ss_alb_bc2(3), ss_alb_bc2(4), &
          ss_alb_bc2(5)
        if (DO_SNO_OC) then
          write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                  &hydrophillic OC: ', &
            ss_alb_oc1(1), ss_alb_oc1(2), ss_alb_oc1(3), ss_alb_oc1(4), &
            ss_alb_oc1(5)
          write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                  &hydrophobic OC: ', &
            ss_alb_oc2(1), ss_alb_oc2(2), ss_alb_oc2(3), ss_alb_oc2(4), &
            ss_alb_oc2(5)
        end if
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &dust species 1: ', &
          ss_alb_dst1(1), ss_alb_dst1(2), ss_alb_dst1(3), ss_alb_dst1(4), &
          ss_alb_dst1(5)
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &dust species 2: ', &
          ss_alb_dst2(1), ss_alb_dst2(2), ss_alb_dst2(3), ss_alb_dst2(4), &
          ss_alb_dst2(5)
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &dust species 3: ', &
          ss_alb_dst3(1), ss_alb_dst3(2), ss_alb_dst3(3), ss_alb_dst3(4), &
          ss_alb_dst3(5)
        write (stdout,*) 'SNICAR: Mie single scatter albedos for &
                &dust species 4: ', &
          ss_alb_dst4(1), ss_alb_dst4(2), ss_alb_dst4(3), ss_alb_dst4(4), &
          ss_alb_dst4(5)
      end if
    end if
  end subroutine SnowOptics_init

  subroutine SnowAge_init( )
   use mod_clm_varctl     , only : fsnowaging
   type(clm_filetype) :: ncid   ! netCDF file id
   ! Open snow aging (effective radius evolution) file:
   allocate(snowage_tau(idx_rhos_max,idx_Tgrd_max,idx_T_max))
   allocate(snowage_kappa(idx_rhos_max,idx_Tgrd_max,idx_T_max))
   allocate(snowage_drdt0(idx_rhos_max,idx_Tgrd_max,idx_T_max))

   if (myid == italk) &
     write(stdout,*) 'Attempting to read snow aging parameters .....'
   call clm_openfile(fsnowaging,ncid)

    ! snow aging parameters

   call clm_readvar(ncid,'tau',snowage_tau)
   call clm_readvar(ncid,'kappa',snowage_kappa)
   call clm_readvar(ncid,'drdsdt0',snowage_drdt0)

   call clm_closefile(ncid)

   if (myid == italk) then
     write(stdout,*) 'Successfully read snow aging properties'
     if ( debug_level > 3 ) then
       ! print some diagnostics:
       write (stdout,*) 'SNICAR: snowage tau for T=263K, &
               &dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_tau(3,11,9)
       write (stdout,*) 'SNICAR: snowage kappa for T=263K, &
               &dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_kappa(3,11,9)
       write (stdout,*) 'SNICAR: snowage dr/dt_0 for T=263K, &
               &dTdz = 100 K/m, rhos = 150 kg/m3: ', snowage_drdt0(3,11,9)
     end if
   end if
  end subroutine SnowAge_init

end module mod_clm_snicar
! vim: tabstop=8 expandtab shiftwidth=2 softtabstop=2
